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ABSTRACT 

We measure the local galaxy far-infrared (FIR) 60-to-lOO |J.m colour-luminosity 
distribution using an all-sky IRAS survey. This distribution is an important reference 
for the next generation of FIR-submillimetre surveys that have and will conduct deep 
extra-galactic surveys at 250-500 iJ.m. With the peak in dust-obscured star-forming ac- 
tivity leading to present-day giant ellipticals now believed to occur in sub-mm galaxies 
near z ~ 2.5, these new FIR-submillimetre surveys will directly sample the SEDs of 
these distant objects at rest-frame FIR wavelengths similar to those at which local 
galaxies were observed by IRAS. We have taken care to correct for temperature bias 
and evolution effects in our IRAS 60 |j.m-selected sample. We verify that our colour- 
luminosity distribution is consistent with measurements of the local FIR luminosity 
function, before applying it to the higher-redshift Universe. We compare our colour- 
luminosity correlation with recent dust-temperature measurements of sub-mm galaxies 
and find evidence for pure luminosity evolution of the form (1 -I- z)^. This distribu- 
tion will be useful for the development of evolutionary models for BLAST and SPIRE 
surveys as it provides a statistical distribution of rest-frame dust temperatures for 
galaxies as a function of luminosity. 
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1 INTRODUCTION 

Deep extra-galactic surveys at 
200-1200 trin) over the last 10 
population of luminous infrared galaxies 



sub-mm wavelengths (~ 
years have uncovered a 
(L > 10^2 Lq) 



with star-formation rat e s inferred to be ^ I OOOMq /yr 
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after) are believed to be high-redshift {z > 1) ana- 
logues, and in many cases more luminous examples, of lo- 
cal Ultra-Luminous Infr ared Galaxies (ULIRGs) discovered 
with IRAS 20 years ago ((Sanders & MirabelllTggg ). Further- 
more, the rest-frame wavelengths sampled by sub-mm sur- 
veys of the highest-redshift SMGs approaches those of the 
far-infrared (FIR) IRAS observations. Appearing in vast 
quantities consistent with massive evolution of the local 
ULIRG population, these SMGs are now believed to repre- 



sent an important early stage in the evolutionary sequence 
that ultimate ly prod uces locally-observed massive elliptical 
galaxies (e.g. IScott et al. 2002; Blain et al., ,2004 ). Thus, to 
this day, the local far-infrared luminosity function measured 
by IRAS continues to be useful for interpreting the results 
of these longer- wavelength surveys. 

Additional motivation for studying the IRAS luminos- 
ity function, and its connection to the higher-redshift SMG 
population, comes from the shape and magnitude of the cos- 
mic infrared backgroun d (GIB) measured by GOBE which 
peaks near 200 \xm. (e.g. lFixsen et al" ] | 19981 ). Its broad shape 
resembles the superposition of many thermal SEDs, which 
can be interpreted as evidence for a population of sources at 
redshifts z < 1 with a large range of physical temperatures, 
or alternatively, as a population with a narrower range of 
temperatures, but residing over a greater range of redshifts, 
including a significant fraction at z > 1 (the SMG popu- 
lation). This latter possibility is su pported by the fact that 
the total energy density of the GIB l|Franceschini et al]|200ll ) 
exceeds the contribution of local I RAS galaxies by a factor 
of ~ 3 jSoifer fc NeugebaueillT99ll) . 
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In this paper we examine the colour-luminosity corre- 
lation and luminosity function of IRAS galaxies, which to- 
gether are an important reference for constraining models 
of galaxy evolution with a new generation of sub-mm sur- 
veys at shorter wavelengths. It has been known for some 
time that more lumino us IRAS galaxies exhibit w armer 
dust temperatures (e.g. ISoifer fc Neugebauej Il99ll ). This 
relationship is important for a class of phenomenological 
models (e.g . iBlain fc Longaiij IT993.: Guiderdoni et al]|l997l: 
Blain et all Il999l: ICharv fc Elbaj I2001I: iMalkan fc Steckeij 



2OOII : iRowan-Robinsonj I2OOII : iLagache et all |2003| . l2o6i ) 
that have been used to predict the source counts 
and redshift distribution s at FIR waveleng ths for the 
Spitzer Space Te l escope jWerner et al.l bOO^ ) and Akari 
ijMatsuhara et aU I2OO6I ). and at sub-mm wavelengths , 
for instru ments such a s SCU BA (|Holland et al.l 19991). 
MAM BO jKrevsa et al.l l2002h. LABOCA (iKrevsa et al.l 
l2003l). gersc/te//SPI RE (iGrifHn et aP |2006|) . SCUBA-2 



( Holland et all l2006h. AzTE C (| Wilson et al.l |2008| ). and 



BLAST (jPascale et all 120081 ). These models often use the 
shape of the CIB as an integral constraint, since the to- 
tal surface brightness of galaxies at each wavelength can- 
not exceed the measured diffuse background. These authors 
apply evolution to local luminosity functions to obtain esti- 
mates of the redshift-dependent luminosity function ^{L, z). 
In order to compare these models with observations, spec- 
tral energy distribution (SED) templates are adopted to ex- 
trapolate observed flux densities from the rest-frame lumi- 
nosities. In a number of cases, it has been beneficial to fit 
data over a range of wavelengths by dividing the local lumi- 
nosity function into several discrete populations, each with 
diffe rent SED templates and separate evolutionary forms 



fe.g.lBlain et al.lll999l :lR owan-Robinsonl200ll : lLagache et al.l 



12003! . I2OO4I ). Recentlv IWall et al.l l|2008h demonstrated di- 
rect evidence for the presence of at least two significant 
populations in a sample of sub-mm luminous sources in 
GOODS-N. The luminosity- colour correlation can be useful 
for such models as a method for assigning dust tempera- 
tures to SED templates as a function of luminosity. For ex- 
ample, [Eagachc ct al. ( 200 3i) use the observed IRAS c olour- 
luminosity distribu tion of lSoifer fc Neugebaueij (|l99ll ), and 
iLewis et al.l (|2005|) use the F IR colour-luminosity distribu- 
tion of Chapman et al.l l|2003l . henceforth C03) inferred from 
a model fit to IRAS data. 

C03 calculate $(L, C), the galaxy volume density as 
a function of total 3-1100 \im luminosity, Lt, and the 60- 
100 n^m colour, C = log(S6o/S'ioo). They formulate <I>(L,C) 
as the product of a luminosity function, and the distribution 
in C as a function of L. The two functions are fit indepen- 
dently, with the latter being constrained directly from the 
observed distribution of the ratio of broad-band IRAS 60 
and 100 vim fluxes. 

In this work we provide a more accurate measurement 
of the joint colour-luminosity distribution using a single 
maximum-likelihood optimization to solve for all of the 
model parameters simultaneously. Our methodology also 
differs from that of C03 in several other key respects: (i) 
rather than calculating C with observed IRAS broad-band 
fluxes we use rest-frame monochromatic 60 and 100 |J.m flux 
densities derived from fltted SEDs; (ii) we use narrower- 
bandwidth 42.5-122.5 [im FIR luminosities instead of 3- 
1100 nm Total Infrared (TIR) luminosities to minimize the 



dependence of the fltted distribution on the choice of SED 
templates; (iii) a correction for redshift evolution in the 
IRAS galaxy population is applied; and (iv) we account 
for a bias against the detection of cooler sources caused 
by the 60 [im selection criterion for the sample. Our galaxy 
sample, SED fltting procedure, and methods for calculating 
luminosities and volumes are described in Section (2] The 
luminosity function and colour- luminosity distribution are 
calculated in Section O We discuss the choice of luminos- 
ity variable, and its consequences, in Section 14.11 Finally, 
we compare the local colour-luminosity correlation with the 
observed values for high-redshift sub-mm galaxies in Sec- 
tion |4]2] and test a simple evolutionary model. Throughout 
this paper a standard cosmology is adopted with Qm ~ 0.23, 
!:!a = 0.77 and Ho = 74 kms"^ Mpc 



2 SAMPLE PREPARATION 

We use the same fl ux-limited Seo > 1.2 Jy IRAS sample of 
iFisher et al.l l| 19951 ) as in C03. Their catalogue covers most 
of the sky and provides 60 and 100 ^.m flux-densities as well 
as spectroscopic redshifts for each galaxy. The cool, high- 
luminosity region of the observed colour-luminosity plane 
found to contain a large number of spurious sources in C03 
has also been excised. We use this sample to flrst calculate a 
non-parametric (binned) FIR luminosity function, and then 
flt it with simple parametric models. 



2.1 SEDs, Luminosities and Colours 

To calculate rest-frame luminosities and colours from ob- 
served IRA S 60 and 100 yu c n flux densities, we follow the 
method of ISaunders etHI (|l990l ). A single temperature 
modifled blackbody SED is assumed for each source, Siv) = 
Au'^ B^{Tobs)- The dust emissivity index is flxed at /3 = 1.5 
which is consistent with typical values measured for local 
ULIRGs with su b -mm follow-up ( e.g. iDunne et al.l |200G| ; 
iKlaas et all I2001I ; lYang fc Phillipg |2007| ). AU of the sub- 
sequent analysis in this paper has also been repeated us- 
ing values 13 — 1.0 and 2.0, and the variation in the re- 
sults is well within the quoted uncertainties. The remaining 
two parameters, the amplitude A, and the observed tem- 
perature Tobs, are then uniquely determined from the ob- 
served Seo and Sioo- Fo r this flt we take into account the 
broad IRAS passbands (jBeichman et al.lll988l ). Bolometric 
luminosities are calculated by integrating the fltted SED di- 
rectly — the bolometric flux emitted in the rest-frame, Shoi 
is simply the integral of the observed SED across the red- 



shifted band, Shoi = f'^,^"r-, , \ S{v)du, where A; = 122.5 vim 
and Au= 42.5 vim for FIR fluxes, and A; = 1100 vim and 
All = 3 vim for TIR fluxes. Similarly the colour C is calcu- 
lated from the logarithm of the ratio of monochromatic flux 
densities emitted at 60 and 100 vim in the rest-frame by the 
model SED. 

Rather than a simple modified bla ckbody, C03 adop t 
the range of model SED templates from iDale et al. I (|200lh . 
This difference has a negligible effect on the inferred FIR 
luminosities sinc e there is very little structure in the 
iDale et all l|200ll ) SEDs at 42.5-122.5 ^m that is not charac- 
terized by the single temperature variable in our SED model 
(despite the correlation between luminosity and P assumed 
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Figure 1. Observed IRAS FIR colour distribution (C = 
log(S'6o/'S'ioo)) as a function of 3—1100 \mi TIR luminosity. Stars 
and triangles show the mean and 68% confidence intervals when C 
is calculated from broad-band IRAS fluxes following the prescrip- 
tion in C03. The dotted lines show the mean and l-cr envelope 
of the fitted COS colour-luminosity correlation for reference. The 
thick solid and dot-dashed lines show the mean and 68% con- 
fidence intervals of the distribution when C is calculated with 
monochromatic flux densities emitted at 60 and 100 (am in the 
rest- frame. 

in iDale et al.l [200lh . For example, assuming temperatures 
ranging from 30 to 50 K (and f3 — 1.5), the FIR luminosities 
inferred from our modifie d blackbody templates compared 
with the lDale et al.l (|200ll ) SEDs with the same correspond- 
ing values of C agree to within ~ 5%. Since the difference 
is so small, and for the sake of simplicity, we therefore pro- 
ceed with the modified blackbody SED model to measure 
the FIR colour-luminosity distribution. The TIR luminosity, 
however, cannot be estimated from the modified blackbody 
model as there is significant emission in the mid-infrared 
(MIR) spectrum (~ 3-60 ^.m) that is misse d by the steep 
drop on the Wien side (e.g. lBlain et al.ll2003 '. and discussion 
in Section [4.ip . Our modified blackbody SEDs in this tem- 
perature raBgejmdeij^^red the TIR luminosities obtained 
from the lDale et al] l|200ll ') SEDs by about 30%. We only use 
their templates to calculate TIR luminosities that are con- 
sistent with C03 for the discussion in this section (Figure [1]), 
and Section [211] (Figure [31 . 

Another fundamental difference between this work and 
CQ3 is the definition of FIR colour. Whereas we choose to de- 
fine C in terms of the ratio of rest-frame monochromatic flux 
densities, C03 use observed broad-band /7?^5 fluxes. We be- 
lieve our definition is more useful as a general reference since 
no detailed knowledge of the IRAS passbands is required 
in order to use our colour-luminosity distribution. Further- 
more, we flnd that the strength of the colour-luminosity cor- 
relation is relatively diluted when using broad-band fluxes. 
In Figure [T] we show the distribution of 77?^ 5 galaxy colours 
as a funct i on of TIR luminosity derived from fits of the 
[Dale et all l|200lh SED templates. The stars and triangles 



correspond to the mean and 68% confidence intervals using 
observed broad-band IRAS fluxes, effectively re-producing 
the top panel of Figure 1 in C03. The thick solid line and 
dot-dashed hues show the mean and 68% confidence inter- 
val of the colour distribution using the ratio of monochro- 
matic 60 and 100 (im fiux densities emitted in the rest-frame. 
At luminosities L <, W^" Lq there is little difference in the 
shapes of the distributions. At L ^ 10^'' Lq , however, the 
colour-luminosity correlation is significantly steeper. 

The dotted lines in Figure [T] show the mean, and 1-a 
envelope of the C03 colour-luminosity correlation. We note 
that although the mean of this fitted parametric distribution 
clearly tracks the stars in the plot, the standard deviation of 
the distribution, ac = 0.065, appears to have been under- 
estimated. We find that both the 68% confidence intervals, 
and the standard deviations of C, for each luminosity bin, 
are typically closer to 0.13. 

2.2 Evolution in the Sample 

Since the most distant, luminous objects in IRAS sam- 
ples exhibit the effects of strong luminosity and/or density 
evolution (e.g. Saunders et al. 1990; Kim & Sanders 1998|; 
[Lawrence et al.li 19991') we must account for its effect in our 
measurement of the local luminosity function. Rather than 
fitting for this evolutionary form ourselves, we instead apply 
an explicit correct ion based on the luminosity evolutionary 
form fit by Sa unders et al. the luminosity of each 

galaxy is divided by a factor (1 -I- z)^ corresponding to its 
redshift. We have chosen to apply a luminosity, rather than 
a density evolution correction for consistency with the dis- 
cussion in Section [4.21 

2.3 Accessible Volumes 

The X;(V^max) estimator (|Schmidtlll96i ). with accessible 
volumes Vmax corresponding to the largest redshift at which 
a galaxy would be detected given the survey fiux limit, used 
in C03, is appropriate for the monochromatic 60 [im lumi- 
nosity function. In this case a given object's luminosity, Lgo, 
is a function only of the observed flux density and distance. 
Therefore the maximum volume in which the given object 
can be detected corresponds to the distance at which the 
observed Sao, given its Leo, drops below the flux limit of 
the sample. For the broad-band FIR luminosity function 
described here, however, the relationship between Lp and 
5*60 is more complicated and this simple method is invalid. 
There exists a bias against the dete ction of cooler source s 
given the shape of the galaxy SEDs (jSaunders et al.lll99Gl ). 
The wavelength of the peak FIR emission is typically in the 
range 60-200 [im. The SEDs of warmer objects peak closer 
to 60 ^tm, and colder objects at longer wavelengths, so that 
in general for a fixed 60 \im flux density a colder object must 
be more FIR lumin o us to be included in the sample. 

ISaunders et~aLl (|199G| ) derive the FIR luminosity func- 
tion from their 60 \im flux-limited survey by selecting a 
sub-sample of objects brighter than a FIR flux limit that 
corresponds to their 60 jim flux-limit and the coolest dust 
temperature that th ey observed, 23 K (see Section 6.5 in 
ISaunders et al.lll990l ). However, this selection reduces the 
size of their sample from ~ 3000 objects, to 1004. There is 
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also an underlying assumption that there is no significant 
population of sources with dust temperatures T < 23K. 

In this work we use the entire sample, but calculate ac- 
cessible volumes using a modified formalism that accounts 
for the dependence of Lp on the FIR colour. Given an ob- 
served temperature, Toba, and redshift, z, the rest-frame 
temperature, T, and total luminosity for an object are calcu- 
lated. The accessible volume corresponds to the maximum 
redshift at which an object with its rest-frame luminosity 
and temperature would be detected in the sample, or corre- 
spondingly the distance at which its observed flux density 
in the IRAS 60vim passband is 1.2 Jy. 



3 THE $(L,C) DISTRIBUTION 

3.1 Non-Parametric (binned) Estimate 

With luminosities and accessible volumes for all of the ob- 
jects in the sample, we first calculate the non-parametric 
(binned) colour- luminosity distribution, $b{L,C). Since the 
accessible volume is now parameterized by both L and C, 
the modified ^(1/Vmax) estimator is simply , 



47r ^ — > 

$b{L,C)dLdC ^ 



where ^^{L, C)dLdC is the number of sources in the area 
of the L-C plane, and the sum runs over all of the galax- 
ies, i, with luminosity-evolution corrected luminosities (Sec- 
tion 12. 2p and colours that land within the bin, and with 
accessible volumes Vi. The factor in front of the sum is the 
fraction of the sky covered by the survey. The binned lu- 
minosity function may be derived from this distribution by 
marginalizing over C, 



<i>biL)dL = $b(L, C)dLdC, 



(2) 



where j runs over all of the bins along the C axis. We define 
a second representation of the luminosity function using a 
lower-case (A, 



0(L) = ln(10)L$(L), 



(3) 



which changes the units from Mpc"'' Lg^ to the more typical 
Mpc~^ dex~^, in order to assist comparison with other work. 
This representation of our FIR luminosity function is shown 
in Figure [2] 

At luminosities L ^ 10® Lq there is excellent agree- 
m ent between our luminosity function and the measurement 
of ISaunders et al] (|l990l ) (shown in Figure [2] with a dotted 
line). However, at fainter luminosities our luminosity func- 
tion includes many more objects. The reason for this discrep- 
ancy is probably due to our choice of ^(l/Vmax) estimator, 
and the fact that over this luminosity range the sample is 
dominated by an over-density of galaxies in the Local S uper- 
cluster. It is for this reason that ISaunders et"al] (|l990D used 
an alternative estimator that is insensitive to local density 
variations. Their method has the potential to more accu- 
rately determine the shape of the luminosity function, how- 
ever at the expense of losing the absolute normalization. 
At luminosities L ^ 10^ Lq the galaxies are typically suf- 
ficiently distant that this issue is no longer important, and 
both estimators give consistent answers (see Section 8 of 




Figure 2. The non-paramotric 42.5—122.5 |am FIR luminosity 
function (symbols with 68% Poisson error bars, and arrows show 
95% upper-limits for bins with < 2 objects) and two parametric 
fits (Equation [6] as a solid line, and Equatio n [7] as a d ashed line) 
derived from the Seo > l-2Jy Jfiyl5 sample oflFisher et al. ( 1992). 
The parametric FIR luminosity function of lSaunders et al.l |199y) 
is shown for comparison as a dotted line. 
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Figure 3. The non-parametric 3—1000 |am TIR, luminosity func- 
tion (symbols with 68% Poisson error bars, and arrows show 95% 
upper-limits for bins with < 2 objects) and the parametric fit of 
COS (solid line). The discrepancies at L > 10^ Lq are due to the 
corrections described in Sections 12.21 and 12.31 



[Lawrence et al.lll999l ). In addition to this effect, lYun et al.l 
(l20(3ll) suggest that some o f the flattening at faint luminosi- 
ties in the Saunde rs et al.l ([T990) luminosity function may 
be due to sample incompleteness. 

We also produce the non-parametric TIR luminosity 
function for comparison with the parametric form of COS 
in Fig ure [S] For this calculation we have used the same 
iDale et al. (2001) SED templates as C03 to derive TIR lumi- 
nosities. The COS modeQ has a significantly different shape 



1 We take the luminosity function to be (Tcj(27r)i''^<I>i(L) from 
Section 3.2 in COS. The dimensionless pre-factor is needed since 
their colour distribution, #2(C), is an un- normalized Gaussian 
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compared with our binned representation at luminosities 
L > lO'' 1/0 , the range over which the local over-density of 
galaxies is irrelevant. Their model under-predicts the binned 
luminosity function by a factor ~ 30% at 5 x 10^" Lq , and 
rises to over-predict by a similar factor at 2 x lO^'^ Lq . We 
have determined that this discrepancy can be explained en- 
tirely by the effects described in Sections 12.21 and 12.31 since 
the model and the bins are otherwise consistent without 
them. Applying only the correction for evolution in the sam- 
ple we find that the number of objects in the brightest 
bins decreases - explaining the factor of 30% at luminosities 
> IQi^ Lq. The reason for this is that these objects are the 
most distant, and therefore exhibit the strongest effects of 
redshift evolution. This correction has almost no effect by 
10^^ 1/0. In contrast, applying the correction for accessible 
volumes increases the numbers of objects in bins at luminosi- 
ties primarily ^ 10^^ Lq. This increase is caused by the fact 
that less luminous objects are cooler, with correspondingly 
smaller volumes in which they could be detected given the 
60 ^-m flux limit. Together, these corrections demonstrate 
that the luminosity function is in fact significantly steeper 
than the result of C03 at luminosities 10^*^ Lq , the most 
important range for comparison with results from new sub- 
mm surveys of distant star-forming galaxies. 



3.2 Maximum Likelihood Model Fits 

Next we fit simple parametric models to the data by max- 
imizing the likelihood of observing the sample. For the re- 
mainder of the paper we consider only FIR luminosities to 
avoid dependence on assumptions about the shape of the 
mid-infrared SED (Section 14. ip . At a given position in the 
L-C plane, the expected number of sources from our sample 
to have landed in that bin, given a model for ^{L, C), is 



fi{L,C) = K.ax?^$(L,C)dLdC. 
47r 



(4) 



These expectations are used to calculate the joint Poisson 
likelihood of the data. 

We express the model, $(L,C), as the product of the 
luminosity function, and the conditional probability 

of a galaxy having a color C given the luminosity L, p{C\L), 

$(L,C) =$(L)p(qL). (5) 



3.2.1 Parametric forms of ^{L) 

For ^(L), we consider two forms. The first is the dual power- 
law of C03, 



(6) 



where L« is the characteristic knee luminosity, p, is 
the number density normalization of the function at L, 
(Mpc~^Lg^), and a and fJ characterize the power-laws at 
the faint {L < L,) and bright {L > L,) ends respectively of 
the luminosity function. 

The other form consider ed is the hybrid p ower- 
law/Gaussian form preferred bv lSaunders et aU l|l99Cf ), 



with standard deviation crc. Also, we have assumed that the units 
for pt are Mpc""^ dex~^ rather than Mpc~''Lg as indicated. 



(l-c) 



exp 



ln(10)I/' 



(7) 



Maximum likelihood solutions for both forms are shown 
in Figure [2] At the faint end (Lfir < L*) both functions 
approach power-laws, and are therefore indistinguishable. 
At the extreme bright-end they diverge; (ps curves below 
at L ~ 10^^ Lq , although both forms lie mostly within the 
error bars of the non-parametric estimate. 

To characterize the quality of the fits we calculate val- 
ues of reduced for luminosity bins that contain at least 10 
objects, approximately luminosities 5 x 10^-10^^ Lq (with 
this number of objects the Poisson error distribution is rea- 
sonably approximated by a Gaussian). Over this range the 
power-law form produces a value of reduced = 2.2, and 
the hybrid form 2.4. Given the similarity of these values, 
and the fact that each form has the same number of pa- 
rameters, we feel there is no compelling evidence to favour 
one model over the other given the data. While the choice 
has no impact on the subsequent discussion in this paper, 
we note that the two forms rapidly diverge at luminosities 
> 10^^ Lq, potentially the most important region of the lu- 
minosity function for comparison with the results of sub-mm 
surveys. For example, while at 10 ^^ Lq the power-law only 
exceeds the ISaunders et al.l l|l990l ) form by about 10%, at 
10^'^ Lq , it is nearly an order-of-magnitude larger. Fitted 
parameters for both models are given in Appendix]^ and 
they should only be considered valid to a maximum lumi- 
nosity of ~ 2 X 10^^ Lq. 



3.2.2 Parametric form of p{C\L) 

In C03 it was shown that the distribution in C is approx- 
imately Gaussian at a particular value of L. The precise 
functional form we have adopted is 



p{C\L) = 



1 



exp 



2 \ ac 

with the mean colour at a given luminosity given bjjf] 



C.-51ogio (i + T;) +7logio (1 + 77) 



(8) 



(9) 



Note that unlike C03, the "knee" luminosity, L' , for p{C\L) 
is independent of the knee luminosity for the luminosity 
function, L*. In addition, the width of the colour distri- 
bution, CTc, is characterized by two different values: af and 
(76 at the faint and bright ends respectively, with a smooth 
transition at L', 



.-L'/L 



-L/L's 



(10) 



Figure |4] compares the mean and l-cr envelope of the para- 
metric p{C\L) with a non-parametric estimate created by 
factoring the smooth model ^{L) from the binned ^hiL, C). 



^ We note that in Section 3.2 of COS the expression for Co is 
clearly meant to be the logarithm of the third equation in Sec- 
tion 3.1 — the form we have adopted here. 
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Figure 4. Comparison between parametric (solid and dotted 
lines give the mean and 1-cr envelope of Equation |8] respec- 
tively) and non-parametric estimates of p{C\L) (greyscale shows 
^h{L,C) normalized along the C axis). The temperature axis is 
derived from C assuming a dust emissivity index /3 = 1.5. At 
luminosities L ^ 10* and L ^ lO^^^ the sample does not contain 
enough galaxies to accurately constrain the shape, and p{C, L) is 
simply extrapolated in the parametric model. 



It is argued in C03 that the width of the distribution in 
C is constant as a function of luminosity. The top panel of 
Figure 2 in C03 demonstrates a constant width in Sm/S\oo 
with a logarithmic axis, in contrast with the bottom panel 
in which the width of the distribution is shown to broaden 
at greater luminosities when plotted with a linear axis. This 
behaviour motivates the definition C = log(S'6o/Sioo)- How- 
ever, this plot appears to be at odds with the top panel 
of Figure 4 in COS which exhibits a systematic broadening 
at higher luminosities. Such a trend does not appear to be 
present in our measurement of p{C\L). For clarity, we com- 
pare slices of the parametric estimate of p{C\L) at several 
fixed luminosities with the binned estimate in Figure [S] The 
good agreement between these two estimates, both in terms 
of scatter and systematic variations, indicates that Equa- 
tion |8] adequately describes the shape of p{C\L). We find 
that the width of the distribution narrows with increasing 
brightness to at = 0.13 from af = 0.2 at a transition lu- 
minosity of ~ 3.5 X 10^ (Appendix IX}. The broadening 
shown in C03 at greater luminosities does not appear to be 
caused by their choice of SEDs or choice of broad-band over 
monochromatic colours. The most likely explanation is an 
artifact of the C03 griding scheme. They use variable-width 
luminosity bins which contain equal numbers of objects, in 
contrast to our method which uses equally-spaced logarith- 
mic bins. The wide, sparsely-populated high-luminosity bins 
may simply dilute the colour-luminosity correlation. We note 
that our fitted value for at appears to be consistent with the 
width of the distribution in the lowest, narrowest, luminos- 
ity bins in the top panel of Figure 4 from C03, (although 
they claim a smaller standard deviation of 0.065; see Fig- 
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Figure 5. The distribution of FIR colours, C, about the mean, 
Co (Equation |9]l, at a range of luminosities. The solid lines 
are normalized slices of the measured (non-parametric) •p(C\V) 
(<l>b(ii C)I^\^{L) - See Equations [T] and O evaluated at C — Co, 
and numbers indicate logjo (Lfir). The parametric model (Equa- 
tion |8} is shown as a dotted line. 

ure [T] and the discussion at the end of Section 12.11 in this 
paper). 

3.2.3 Parameter Uncertainties 

We characterize the uncertainties in the ten parameters of 
$(L,C), for both parametric forms of $(L), using a boot- 
strap Monte Carlo technique. First, 100 realizations of the 
1.2 Jy survey are created from the actual survey data by 
randomly sampling sour ces from th e catalogue with replace- 
ment (see Section 6.6 of I Wall fc Jen kins 200j). We then fit 
the model to each simulated sample. From these 100 fits we 
calculate the sample variances, and covariances between all 
pairs of parameters to estimate the full parameter covariance 
matrix. The maximum likelihood values of each parameter, 
their standard deviations, and the Pearson correlation ma- 
trices are given in Appendix [X] Note that the parameters 
for p(C\V) are largely independent of the parametric form 
chosen for 3>(L). 



4 DISCUSSION 

4.1 Choice of Bolometric Luminosity Variable 

In this work we have chosen to use the 42.5-122.5 vim FIR 
luminosity instead of the wider-bandwidth 3-1100 ^im total 
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Wavelength [^tm] 

Figure 6. The peak-normalized spectral energy distribution of 
Arp 220 (dotted line), and shifted to z = 2.5 (solid line). The 
complete observed BLAST and SPIRE band (200-600 ^mi, in- 
cluding the 30% finite bandpasses for each channel) is indicated 
by the solid grey shaded rectangle. The rest-frame 3-1100 \xm to- 
tal infrared (TIR) flux, and 42.5-122.5 \xm FIR flux correspond 
to the integrals of the horizontal-line filled, and cross-hatch filled 
regions respectively, for the z = 2.5 SED. Clearly for objects at 
this redshift, the BLAST and SPIRE filters more closely match 
the rest-frame FIR fiux than the TIR flux. The peak-normalized 
M82 SED redshifted to z = 2.5 (shown as a dashed line) illus- 
trates the relatively large variations in the rest-frame mid-infrared 
spectrum (~ 3—60 \im observed wavelength). Symbols indicate 
the wavelengths accessible from ground-based sub-mm surveys 
(e.g. SCUBA at 850 [jm, and BOLOCAM/AzTEC/MAMBO at 
1.1mm), and existing space-based FIR data (e.g. Spitzer at 24, 
70 and 160 l^m), demonstrating their inability to accurately con- 
strain the bolomctric FIR luminosity for these galaxies. 

infrared luminosity, Lt, as in COS. An argument for using 
the latter is that it includes a significant fraction of the total 
power emitted by a galaxy that is missed at shorter wave- 
lengths (< 40 li-m) — an effect which becomes increasingly 
important at high luminosities given the positive L-C corre- 
lation. By using Lpm, the "clipping" of shorter- wavelength 
light obscures the physical interpretation of the correlation 
between the two variables. 

Our primary goal, however, is to assist with the de- 
velopment of a model for the luminosity f unction of SMGs 
discovered in BLAST JPascale et al.ll2008l ) and future Her- 
schel/SVmE (|Grif5n et all |2006| ) 250-500 urn surveys, as 
well as any other surveys at similar wavelengths. It is now 
generally accepted that the redshift distribution for the bulk 
of SMGs discovered in 850 \im SCUBA surveys peaks at red- 
shifts ~ 2 .5. For example, the radio-d etected spectroscopic 
sample of lChapman et all (|2003l . |2005| ) finds a median red- 
shift of 2.2 with an interquartile range z = 1.7-2.8, in gen- 
eral agreement with several other studies using radio-FIR or 
radio and 24 yim guided photometric redshif t estimates (e.g. 
lAretxaga et al. 1 120031 . l2007l : |Pope et al.|[200^ . Since the neg- 
ative /f-correction produces a nearly un-biased detection ef- 
ficiency at 850 \ira for a typical ULIRG SED at 2; ~ 1-8 (e.g. 



iBlain et al.]|2002l ), the observed redshift distribution is a 
reasonable proxy for the total dust-obscured star-formation 
rate history of massive galaxies. If SMGs have thermal SEDs 
similar to the ULIRGs that populate the bright end of the 
IRAS luminosity function presented here, their SEDs peak 
at wavelengths ~ 60-200 \im in the rest-frame, and they are 
redshifted into the 200-600n.m BLAST and SPIRE band- 
passes near the peak of their redshift distribution. 

As an example. Figure [S] shows the SED of the ULIRG 
Arp 220 at a redshift z = 2.5 compared to the BLAST 
and SPIRE bandpass region, and the integrated FIR and 
TIR fluxes. Clearly the BLAST and SPIRE integrated fluxes 
more closely matches the rest-frame FIR than TIR flux. 

The peak-normalized SED of the starburst galaxy M82 
redshifted to z = 2.5 is shown for comparison as a dashed 
line to illustrate the relatively large scatter at mid-infrared 
(MIR) wavelengths compared to the smooth thermal SED at 
longer wavelengths. The rest-frame MIR SEDs of these ex- 
ample galaxies exhibit prominent polycyclic aromatic hydro- 
carbon (PAH) absorption and emission features (the ~ 100- 
10 [im range of the observed SED; measurements for actual 
SCUBA-selected SMGs are given in iPopc ct al. 2008). For 
these two examples there is a difference of ~ 30% in the 
contribution of the mid-infrared emission to the TIR lumi- 
nosity. For these reasons, evolutionary models based on the 
FIR luminosity function are less dependent on assumptions 
about the intrinsic near-IR to millimetre-wavelength SEDs 
of high-redshift galaxies than the TIR luminosity function, 
enabling cleaner comparison with new and future data from 
deep sub-mm cosmological surveys. We emphasize the fact 
that one is free to adopt any template library to infer flux 
densities at other wavelengths for objects drawn from our 
$(Z/, C) distribution provided that they have roughly ther- 
mal FIR spectra, and span the relevant range of rest-frame 
FIR colours —0.65 ^ C ^ 0.25 (see for example the com- 
parison in Section r2.1l be tween our modifie d blackbody SED 
model and the library of lDale et al.ll200ll ). 

4.2 Evidence for an Evolving Colour-Luminosity 
Correlation 

Evolution in the FIR colour-luminosity distribution is dif- 
ficult to probe from IRAS catalogues given their relatively 
low redshifts. For exam ple, the median galaxy redshift in 
the lFisher et al.l (Il995l ) sample is z — 0.019, and the most 
distant object is at 2; = 0.326. Despite this, it has been pos- 
sible to place weak constraint s on the evoluti o nary form of 
the FIR luminosity function. ISaunders et al.l (|l990l ) found 
that the most distant galaxies could undergo extremes of 
pure density evolution of the form (1-1- z)^*^, or pure lumi- 
nosity evolution of the form (1 + z)^^^ (we use this latter 
form explicitly to correct our sample, see Section [2TT}. Since 
the most distant IRAS galaxies are also the brightest in the 
sample (luminosities > L,) it is not possible to determine 
which form (or combination) is the more relevant. 

At higher redshifts, the best constraints on the evolu- 
tion of ultra-luminous infrared galaxies come from SCUBA 
surveys at 850 \im. It can be shown, for example, that con- 
tinued luminosity evolution in the bright-e nd of the local 
FIR lu minosity function of the same form as lSaunders et al] 
1 19901 ) to redshifts ~ 2-3 can be used to fit mm-submm 
source counts, assuming ULIRG-like SEDs to extrapolate 
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Figure 7. Comparison of local p{C\L) (shaded region is the 68% 
confi dence interval , solid line is Co from Equation ^ with data 
from lCoppin et al. 1 ||2008|) . The temperature axis is derived from 
C assuming a dust emissivity index /3 = 1.5. Stars with dotted 
l-CT error bars indicate the 10 SMGs with spectroscopic redshifts 
and temperatures derived from observed 350 |am/850 \xm colours. 
Squares have been drawn around the symbols for the 6 objects 
at 2: > 2. High-redshift ultra-luminous galaxies appear system- 
atically cooler than those in the local universe. Under the as- 
sumption of pure luminosity evolution of the form (1 -|- 2)^, the 
SMGs have been projected into the local colour-luminosity dis- 
tribution by shifting them along the luminosity axis (diamonds 
with solid error bars). The model is a plausible fit to the data ex- 
cept for LOCK 850.4 and LOCK 850.41 which appear much cooler 
(shown as lighter symbols). These objects may have ambiguous or 
incorrect optical/IR counterpart identifications (discussed in Sec- 
tion gjj. Future BLAST and SPIRE surveys will constrain FIR 
luminosities of ~ 10^^ Lq galaxies to ~ 20%, with uncertainties 
in C of ~ 0.1 (several Kelvin). A representative measurement in 
the L-C plane is given by the thick cross-filled square. 



observed flux densities from r est-frame luminosities (e.g 
IScott et al]|2002l : lLagache et al.| [2003). A more direct test of 
evolution for a small sampl e of SMGs with k nown redshifts 
was recently performed by I Wall et al.l (|2008l ) finding simi- 
lar results. Unfortunately, until recently, shorter-wavelength 
data that would help to constrain the dust temperatures of 
objects in these samples is generally unavailable, and it is 
therefore not possible to search directly for colour evolution. 
Measuring dust temperatures for large samples (S> 1000) of 
SMGs is one of the primary science goals of BLAST and 
Herschel/BVYREi surveys — with the caveat that redshifts 
must first be determined independently for at least a subset 
of these new objects since there is a potential degeneracy be- 
tween the appare nt observed temperature and redshift (e.g. 
iBlain et al.l(2003l ). 

Recent SHARC-II 350 \im o bservations of SMGs (e.g. 
iKovacs et al.|[200^ : ICoppin et al . 2008), however, have en- 
abled improved estimates of dust temperatures for smaller 
samples (several tens of galaxies) . In Figure [7] ten objects 
(triangles) with constrained dust temp eratures and spectro - 
scopic redshifts (median z = 2.1) from ICoppin et al.l (l2008l) 



are com pared with our loca l me asurement of p(C|L ). As 
noted in iKovacs et all (120061 ) and lCoppin et al.l l|2008l ). the 
rest-frame temperatures for such luminous galaxies (median 
Lfir = 2.3 X 10^^ Lq) are much lower than objects in the 
local Universe. With the correlation between luminosity and 
FIR colour in-hand, we now ask the question: can pure lu- 
minosity evolution of the form (1 + z)^ account for the ap- 
parently cooler temperatures of SMGs at high-redshift? Ex- 
plicitly, we express the redshift evolution of the FIR colour- 
luminosity distribution, "1>(L,C, z), as a simple function of 
the local distribution. 



$(L,C,2) = $ 



L 



(1 + ^) 



(11) 



For comparison, we project each observed SMG into the 
local p{C\L) distribution by dividing their luminosities by 
(1 + z)"', shifting the objects to the left in Figure [7] (dia- 
monds). Given the uncertainties, the 8 warmer objects (top 
of the plot) are roughly consistent with the local distribu- 
tion once we apply this transformation. Ignoring the red- 
shift uncertainties and adding the measured colour uncer- 
tainties, (Tdata, in quadrature with the intrinsic colour width 
in the source population, ctc (Equation 1 101 at the evolution- 
corrected luminosity of the galaxy - the shaded region of 
the figure), we calculate residuals, R, between the model, 
Co (solid line through the centre of the shaded region) and 
the measured colours, Cdata: 

Cdata ~ Co 



R = 



(12) 



We note that for five of these objects R < 1, and the re- 
maining three are in the range R = 1-1.5; the approximate 
expectation for uncorrelated Gaussian uncertainties. This 
result is contrary to the conclusion of C03 who found no 
need to invoke redshift evolution for p{C\L) when compar- 
ing with lower-redshift (z < 1) samples. 

The two most significant outliers, LOCK 850.4 (C = - 
2.2, z = 0.526) and LOCK 850.41 (C = -2.47, z = 0.689) , also 
the two coolest and least luminous objects, deserve further 
explanation. These two galaxies are also the only objects 
with redshifts 2 < 1. Since their observed sub-mm colours 
are otherwise similar to the other galaxies, these low red- 
shifts also imply lower rest-frame dust temperatures (lower 
values of C) . There is a possibility that the true optical coun- 
terpart (and hence redshift) for LOC K 850.4 is at z = 1.482 , 
rather than 0.526 as adopted by ICoppin et al.l ||2008|). 
Both potential counterparts were proposed in Ivison et al.l 
(2005). Adopting the higher redshift object as the coun- 
terpart, the inferred FIR luminosity increases from 8 x 
10^" Lq to 10^^ Lq and the rest-frame dust temperature 
from 13 K to 21 K, or C = -l, at which point it would ap- 
pear to have similar dust properties to the other galax- 
ies in the sample. There is a similar possibility of a mis- 
identifica tion for LOCK 850 .41. Two counterparts are sug- 
gested in llvison et al.l (|2005l ) , although they were only able 
to obtain the spectroscopic redshift indicated above for one 
of them. The other counterpart has an optical ph o tomet - 
ric redshift estimate of z = 2.2±0.2 from iDve ets^ (|2008l ). 



A similar val ue of z = 2.1±n. e, base d on its FIR colours 
is proposed in lAretxaga et al.) l|2007l ) . Adopting a redshift 
of 2 = 2.2, the luminosity for LOCK 850.41 increases from 
4 X 10^" Lq to 5 X lO'^^ Lq and the rest-frame dust tempera- 



The local FIR Galaxy Colour- Luminosity Distribution 9 



ture from 12 K to 22 K, or C = -0.9, also closer to the distri- 
bution for the other objects in the sample. However, if the 
lower-redshift candidates for these two galaxies are correct, 
and a population of galaxies with temperatures T ^ 15 K do 
exist in abundance at redshifts z < 1, it is possible that they 
were completely missed in IRAS surveys, and will appear in 
the wide-area SCUBA-2 and BLAST and SPIRE surveys. 

Finally, we note that the higher-redshift objects gen- 
erally fall the closest to the colour-luminosity distribution. 
To emphasize this fact we draw squares around the 6 ob- 
jects in this small sample at redshifts z > 2. Naively this 
fact is slightly surprising, since one might suppose that 
the nearer objects are in fact more similar to the sample 
used to constrain the local distribution. However, an addi- 
tional consideration is the selection function for this SCUBA 
sample. While the negative if-correction produces approx- 
imately the same observed 850 \im flux density at redshifts 
z ~ 1-8 for a fixed FIR luminosity and temperature, there 
is also a bias towards the detection of cooler objects at a 
fixed redshift and fiux density since such objects are less 
luminous, a nd hence more abundant in the rest-frame (e.g. 
lEales et al.l lT999: Blain ct al. 2002; Cha pman et al. 2003). 
For the lower-redshift objects, at which point the negative 
Jf-correction is diminished, the luminosities are also fainter 
for a given flux density, and this bias could be increased due 
to the broadening in the colour-luminosity correlation that 
we have measured at lower luminosities. 

This comparison is by no means an exhaustive study 
of evolution in the FIR colour-luminosity distribution. It is 
our goal to extend this investigation to much larger sam- 
ples of SMGs with FIR colour information in new BLAST 
extra-galactic surveys (Devlin et al. in prep.), and fu- 
ture //ersc/iei/SPIRE surveys. Understanding the details of 
this evolution is intimately related to the star-formation 
rate history of massive galaxies, since the rest-frame FIR 
luminosity is the key observed quantity in SMGs from 



which star-formation rates are derived fe.e.lKenn cutti 


1998| 


iHuEhes et al.lll998l: IScott et al.ll2002l: IPodc et all 2006 


)• 



5 CONCLUSIONS 

We have measured the local FIR galaxy colour-luminosity 
distribution based on the 60 an d 100 vim flux densities and 
redshifts from the iFisher et al ] l|l995.) all-sky IRAS sample. 
This distribution is an important reference for forthcom- 
ing BLAST and //ersc/ieZ/SPIRE extra-galactic surveys at 
250, 350, and 500 [im that will detect thousands of SMGs 
at redshifts z > 1. Since the space density of SMGs appears 
to peak at redshifts z ~ 2.5, deflning the epoch at which 
most of the stars in present-day massive galaxies formed, 
the BLAST and SPIRE bandpasses will sample the same re- 
gion of the rest-frame SEDs for these objects as IRAS at 60 
and 100 \im for local samples. Our measurement is therefore 
the primary present-day boundary condition with which any 
evolutionary model for the luminosity and dust-temperature 
distribution of SMGs must be compared. This applies to cur- 
rent BLAST and future SPIRE surveys, as well as higher- 
resolution ground-based observations with SCUBA-2 for ex- 
ample. 

Our method accounts for a temperature bias in the un- 
derlying 60 fim flux-limited sample, as well as luminosity 



evolution. These corrections indicate that the bright-end of 
the luminosity function is significantly steeper than an ear- 
lier calculation by C03 which neglected them. We have veri- 
fied that our di stribution is consistent with the FIR luminos- 
ity function of ISaunders et al] (|199G| ) by marginalizing over 
colour. We fit a parametric model to the data consisting of 
the product of the luminosity function, ^{L), with the condi- 
tional colour probability, p{C\L), where C = log(S'6o/S'ioo). 
We fit p{C\L) using a normal distribution for C as a func- 
tion of L, with a mean colour given by a broken logarithmic 
function also of L. This is characterized by the knee luminos- 
ity for the break, L' , which is independent of the luminosity 
function knee, L*, and the width of the distribution, by two 
different standard deviations ab and cr/, above and below 
the knee L' . 

We have checked directly for evolution in the colour- 
luminosity correlation using observations of high-redshift 
SMGs (z > 1) with temperatures constrained by 
SCUBA 850 L i m an d SHARC-II 350 ^m photometry from 
ICoppin et al.l l|2008l ). These high- 2; ultra-luminous objects 
appear much cooler than local galaxies of comparable lumi- 
nosities, and there is preliminary evidence that pure lumi- 
nosity evolution in the local colour-luminosity distribution 
of the form (1 + z)^ is consistent with uncertainties in their 
measured redshifts and colours. This result is contrary to 
C03 who find no evidence for a change in the relationship 
between luminosity and colour in low-redshift {z < 1) sam- 
ples. 
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Table Al. Maximum likelihood parameter values, and l-cr un- 
certainties for ^(L, C) (Equation[5]l using the COS dual power-law 
form of the luminosity function (Equation |6l and Equations ISl llOl 
for p(C|L)). 



Parameter 


Value 


P* 


(1.22 ± 0.24) X 10"" Mpc-^Lg^ 


a 


2.59 ± 0.03 




(5.14 ± 0.39) X 10^° Lq 


P 


2.65 ± 0.05 




0.128 ± 0.003 




0.20 ± 0.01 




-0.48 ± 0.02 


5 


-0.06 ± 0.02 


7 


0.21 ± 0.01 


L' 


(3.2 ± 1.7) X 10^ L0 



Table A3. Maximum likelihood paramete r values, and 1-cr uncer- 
tainties for $(L, C) (Equation|5} using the lSaunders et"al ] l|l990l ) 
form of the luminosity function (Equation [3 and Equations ISl llOl 
for p(C|L)). 



Parameter 


Value 


P* 


(6.29 ± 0.64) X 10-^ Mpc-^dex-i 


a 


1.59 ± 0.03 




(3.99 ± 0.53) X lO^ Lq 


a 


0.60 ± 0.01 




0.127 ± 0.004 




0.20 ± 0.01 




-0.47 ± 0.03 


5 


-0.05 ± 0.02 


7 


0.22 ± 0.01 


L' 


(3.8 ± 2.5) X 10^ L0 



APPENDIX A: MODEL PARAMETERS AND 
UNCERTAINTIES 

Maximum likelihood parameters and Icr uncertainties for 
the local FIR colour-luminosity distribution are given in Ta- 
bles |X]l andRti. Pearson correlation coefficients for the un- 
certainties are given in Tables Kb andlXll. 
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Table A2. Parameter Pearson correlation matrix for $(L,C) using the COS dual power-law form of the luminosity function (see 
Table |X]L). 





p* ct 


L* 


/3 






C, 


<5 


7 


L' 


p* 


1.00 -0.88 


-0.93 


-0.24 


0.03 


0.03 


-0.02 


0.04 


0.06 


0.05 


a 


1.00 


0.73 


-0.13 


0.05 


0.00 


-0.04 


-0.06 


-0.18 


-0.13 


L* 




1.00 


0.49 


-0.10 


-0.02 


0.08 


-0.00 


0.04 


0.03 


P 






1.00 


-0.23 


0.03 


0.17 


0.10 


0.31 


0.22 










1.00 


0.09 


-0.77 


-0.74 


-0.71 


-0.76 












1.00 


-0.45 


-0.29 


-0.39 


-0.47 


a 












1.00 


0.92 


0.75 


0.94 


5 














1.00 


0.67 


0.83 


7 
















1.00 


0.89 


L' 


















1.00 



Table A4. Parameter Pearson correlation matrix for <J>(L,C) using the lSaunders et al.l l ll990l) luminosity function (see Table lAfe). 





p* ce 


L* 


a 






C, 


5 


7 


L' 


p* 


1.00 -0.73 


-0.98 


0.83 


0.27 


0.33 


-0.34 


-0.33 


-0.31 


-0.38 


a. 


1.00 


0.75 


-0.34 


0.09 


-0.11 


0.01 


0.07 


-0.09 


0.02 






1.00 


-0.85 


-0.23 


-0.30 


0.29 


0.29 


0.29 


0.35 


a 






1.00 


0.33 


0.31 


-0.33 


-0.29 


-0.40 


-0.41 










1.00 


0.41 


-0.78 


-0.67 


-0.83 


-0.77 












1.00 


-0.74 


-0.67 


-0.55 


-0.67 


c. 












1.00 


0.92 


0.83 


0.93 


5 














1.00 


0.73 


0.81 


7 
















1.00 


0.91 


L' 


















1.00 



